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SUMMARY 


A variational finite element model for transonic small disturbance 
calculations is described. Different strategy is adopted in subsonic and 
supersonic regions, and blending elements are introduced between different 
regions. In the supersonic region, no upstream effect is allowed. If rec- 
tangular elements with linear shape functions are used, the model is similar 
to Murman’ s finite difference operators. Higher order shape functions, non- 
rectangular elements, and discontinuous approximation of shock waves are also 
discussed. 


INTRODUCTION 


The plane^ steady, inviscid flow past a smooth configuration near sonic 
speed can be described by a perturbation velocity potential cf) satisfying 
the transonic small disturbance equation (TSDE) 


(K-c() )(|) + c|) = 0 

X XX yy 


( 1 ) 


where K is a similarity parameter. This equation is nonlinear and of mixed 
hyperbolic-elliptic type. Its weak solution admits discontinuity in the 
pressure, 


<K-(J) > = - 



( 2 ) 


and 



I'"*]) 


I*]] - 0 


(3) i (3’) 


where < > and [[ I signify the average and the jump across the shock x^(y). 
The flow field solution is required to determine the pressure distribution on 
the airfoil (unlike the methods of singularities, or Kernel methods, used for 
incompressible flow calculations). Recently, finite difference solutions 
have been obtained with marked success (refs, 1-4). 


In this paper, the feasibility of applying a finite element approach to 
transonic flow problems will be studied. A finite element method should handle 
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the same problems that finite differences did, namely, the change of the type 
of the equation in the domain of interest with a discontinuous solution 
satisfying prescribed jump conditions. In passing, the potential solution is 
completely reversible (no entropy changes), and an expansion shock must be 
excluded (using an artificial viscosity or a shock fitting procedure) . Hope- 
fully, complicated boundary conditions will be handled easily in the physical 
space, and the use of higher order shape functions will be efficient. 


FINITE ELEMENTS - BACKGROUND 


Elliptic Problems 


(4) 


(5) 


( 6 ) 


If linear shape functions on triangular elements are used, the algebraic 
equations for the nodal values are identical to those obtained by applying a 
centered difference scheme. 

A gradient method for solving this problem is 

= - p 61 (())'') (7) 

where n indicates the iteration and the optimum p may be obtained in terms 
of the Residual and the Hessian. 

Many nonlinear elliptic problems are solved iteratively by casting them 
in Poisson’s form, where nonlinearity acts as a driving force (incompressible 
sources) = 


Consider the classical boundary value problem, 


L (d)) = Kd) + d) = f on a rectangular Q, 
e XX ^yy 


where c() is known on 3fl; K > 0. The associated functional is 

I(4>) = 2fcj) dxdy 

^ ^ 

The first variation is set equal to zero 

1 1 




■ II ^ 


0 0 


(b - f) 6d)dxdy = 0 

yy 


and the second variation is positive definite. 


+ 64> = - ooR(4)’^) (8) 

XX yy ^ 

and where R is the Residual and oj is a relaxation parameter. 

Argyris (ref. 5) calculated compressible subsonic flows by the Galerkin 
method and obtained impressive results within a few iterations. Similar 
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applications were reported by Gelder (ref. 6), Norrie and DeVries (ref. 7) 
Periaux (ref. 8), and Chan and Brashears (ref. 9). 


Hyperbolic Problems 


Finite element methods were also developed for approximate solutions of 
initial value problems. Both variational and weighted Residual methods were 
used (see refs. 10 - 16). Most of these investigators used either finite 
element in space with finite difference in time, a quasi-variational principle, 
or a convolution bilinear form. A variational formulation for initial value 
problems is not possible in the classical context of the calculus of varia- 
tions. Consider the simple linear wave equation 

U (cj)) = K(f) + c|) = f, K < 0 (9) 

h XX yy 

where (j)(x = 0) and ({> (x = 0) are given as initial conditions and x is 
the time-like coordinate. Application of Hamilton’s principle requires know- 
ledge of the conditions at the beginning and end of a time interval and does 
not apply here. This is difficult because we persist in employing boundary 
value techniques to solve an initial value problem. 

Contrary to the conventional shooting method (an initial value technique) , 
which employs a marching (step-by-step) scheme to solve a boundary value 
problem, here we will solve the initial value problem by a formal application 
of Hamilton’s principle. The success of the shooting methods depends on the 
assumption that a variation in the initial slope has a one-to-one correspond- 
ence with a variation in the end position; hence, the problem can be solved 
iteratively. At each iteration, only an initial value problem is solved. For 
linear problems, iterations may not be needed. The reverse of this process is 
valid if the same assumption holds, namely, initial value problems can be 
solved iteratively, with each iteration consisting of a boundary value problem. 
Again, iterations may not be needed for linear problems. 

So, if we assume that the end value c()(x = X) is known instead of the 
initial slope (j) (x = 0), the associated functional (potential and kinetic 
energy) would be^' 

X 1 

= f f = 4'y^ + 2fcj) dxdy (10) 

0 0 

which can be discretized and expressed as a sum over finite elements. A basic 
requirement for application of Hamilton’s principle is that we not vary the 
extreme positions of the physical system. The missing equation (the variation 
with respect to the end position) is replaced by an equation prescribing the 
variation with respect to the initial slope (see fig. 1). 

Note that the second variation is not positive (stationary but not 
extremum) , and there may be no advantage over weighted Residual methods with 
a sensible choice of suitable weighting functions. We note also, that 
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arrangement of the elements is not completely arbitrary, and sometimes 
the element size is restricted by stability requirements. For example, if we 
use linear shape functions on triangular elements, the algebraic equations 
for the nodal values are identical to those obtained by explicit centered 
difference schemes. These requirements arise because a hyperbolic system has 
characteristics (or preferred directions of propagation) and by just minimizing 
the energy, we have not taken these features into account. Implicit 
(unconditionally stable) schemes will be discussed below. 

For many nonlinear hyperbolic equations, the following iterative 
procedure can be used: 

- + 6(|)yy = - WR(<|)") (11) 

where a is determined to guarantee convergence of iterations (the approximate 
domain of dependence contains the exact one). 


TRANSONIC FLOWS 


Consider the functional 


!(<()) 


=JJ' ^ (K4>^ + ‘f’y) - f 'I’^dxdy - f g4> ds (12) 


Perturbing in any direction n(H is an admissible function) 

1(4) + en) = I(4>) f ~ j' 

.2 

I l(K ~ <t> )n^ + 1. 

y 


2 3 /* /* 


Vanishing of the first variation gives 


/ /(K<j> - ^ <|5^)n + 4> n dxdy 

y y X 2 X X ” 


y y 


-/ 


gn ds = 0 


(lA) 


Applying Green's theorem, equation (3) becomes 


//K - i 


(|)^) + ((j) ) ridxdy 

XX y y ^ 


■/ 


(4> - g)n ds 

n 


= 0 (15) 


Note that the second variation is not always positive. 
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Iterative Procedures 


For a nonlinear problem, we need a linearization procedure and a 
discretization technique. In general, they do not commute. 

If we start by discretizing the integral expression, minimization will 
lead to a nonlinear system of algebraic equations to be solved iteratively 
(e.g., Newton’s method). On the other hand, consider the sequence of 
functionals 

At each iteration, only a linear system of equations will be solved. 

Discretization Procedures 


The finite element method has been used to solve efficiently subsonic 
flow problems, with complex geometries employing nonrectangular elements, 
with a better approximation of the boundary conditions than finite differences. 
Although the matrix for the nodal values will not have the same regular 
structure as in finite differences, the number of unknowns is usually less 
(for higher order elements), and the matrix inversion procedure is different 
(banded Gaussian Elimination). 

For transonic small disturbance theory, the streamlines are almost 
parallel to the x-axis, and the body boundary condition can be applied at 
y = 0. Moreover, in the supersonic bubble, x is the time-like coordinate, 
and the nodes may be located along x = constant lines. Finite differences 
suit the problem very well. The small disturbance simplifications eliminate 
the advantages of finite elements. The situation will be different, however, 
if the full potential equation is considered where the flow direction is 
unknown and if the exact boundary conditions are applied at the surface of the 
body. 


Nevertheless, we will consider a simple example and use rectangular 
elements to study the feasibility of using a finite element approach to a 
mixed type equation. As a matter of fact, efficient finite difference schemes 
for elliptic and parabolic equations are constructed this way (see refs. 17 - 
19). 

Semi-Discretization 

Let 

m 

X (x) Y (y) 
i=l ^ ^ 

where m is the number of strips in the y— direction. The functional I 
becomes ^ 
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A A ( fl fl dx. dX. ) 

L ((p) = 2^ 2^ )K.. / X.X. dx + / M, dx> 


where 


K 


y=y. 

/ 


dY. dY. 


IJ dy dy 

y=y^ 


dy 


and 


M 


y=yx 

/ 


(K-4>“) Y.Y.dy = / KpY.Y.dy 
IJ ^ ^ J •'y=0 ^ ^ J 


y=y^ 

/ 


The kinetic and potential energies are 


(17) 


, dX. dX. T 

= iy'y'M . . -A V = .X.X. 

2 ii dx dx 2 ij 1 i 


The Euler-Lagrange equation reads 


6 y*(T-V)dx = 0 (i.e., - (MX) + KX = o) 


(18) 


where M and K are the mass and the stiffness matrices. Or, in the 
canonical form, 


MX = P P = + KX 


(18’) 


where X. (x) must satisfy the essential boundary conditions. For local 

hyperbolic regions, the end value X.(x = x ) will be replaced by an initial 

dX. ^ 

condition, — (x = x. ). 

ax 1 


Full- Discretization 


Instead of solving a system of ordinary differential equations along lines, 
we will consider different discretization procedures also in the x-direction. 

Finite Element in Space, Finite Difference in Time. - If linear hat 
functions in y are used, M will be a triadiagonal matrix 

1^1 4 ij and K will read ^-1 2 J ^ 

These two matrices will be modified by introduction of the boundary conditions. 

In the x-direction, centered differences in the subsonic segment will give 
star A, as shown in figure 2, while backward differences in the supersonic 
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segment will give star B. At the parabolic point P, is set equal to zero. 

At the shock point S the locally normal shock relation <K-cJ) > = 0 provides 
(() downstream of the shock and is used as a derivative boundary condition for 
tSe rest of the unknowns on the line. 


Finite Element in Space and Time . - If linear hat functions in both y 
and X are used, both stars A and B will be the same as in figure 3. 

Higher Order ^h^e Funct ions: Linear Hat Functions in y and Hermite 

Cubic s in x. - The cubic polynomial on 0 ^ x j< Ax , which takes on the four 


prescribed values 

‘f’o’ '^>Xq’ 

(p^, and (|)^^ 

, is 





{X.} = 

1 

o 

o 


01’ 

{ 

=^0 

} 

with 

«00 = 

20^ - 

30^ + 

1 

{ 

4> 

} 


«01 = 

- 20^ 

+ 30^ 


{ 

h 

} 


O 

II 

(6^ - 

II 

(N 

CD 

<N 

0)Ax 

{ 

% 

} 


«ii = 

(0^ - 

0^)Ax 


- 

1 

— 


(6 = 

— ) 







In the subsonic region, the contribution of the neighboring element will be 
included through the assembly of the elemental expression into the global 
system (see figure 4). In the supersonic region, the stationary value with 
respect to cj) (0) and cf) (Ax) ,, assuming 0(0) and 0(Ax) are known, will give 
two algebraic^equations tSat will be used to solve for 0(Ax) and 0 (Ax) 
(according to the inverse shooting method described earlier), namely,^ 


or 




{X} dx = 0 


( 20 ) 


r *11^21 ['*1 T .['ll ®12] [‘♦on^p'on (20.) 

^"^21 ^ 22 -* L ®21 ® 22 J ^*1 


Note, no upstream effect is allowed in the supersonic region, 

Nonrectangular Elements . - All the previous approximates were special 
cases of tensor products. To relax this restriction, consider 



where are the global shape functions. (For example, the isoparametric 

element with four nodes, where <p, and (j) are given at each node, curved 

boundaries are allowed with the restriction thZt the nodes in the supersonic 
region lie on. x = constant lines.) 

Element Equations and Assembly Procedures 


For simplicity, consider a bilinear element with four nodes: 

(j)^ = a + bx + cy + dxy 

The coefficients a, b, c, and d are given in terms of the four nodal values. 

(If the elements were rectangular, this case would reduce to the tensor product 
of linear hat functions in x and y . ) If we consider the element equations 
rather than the nodal equations, the usual finite element assembly procedure 
in the supersonic region must be modified according to the inverse shooting 
method, as shown in figure 5. 

The transition between the elliptic and hyperbolic parts of the flow is 
achieved by introducing blending elements between different regions. Two such 
elements are used: one for the sonic line; one for shock waves. 

Sonic Elements 

For sonic elements, the average of (K-<j)^) is set to zero. These 
elements act as a "buffer zone" between subsonic and supersonic elements. We 
can show that the system matrix will be positive definite if the above 
assembling strategy is adopted and if the sonic element is included. 

Shock Elements 

In transonic small disturbance calculations by finite differences, shocks 
are either captured (using artificial viscosity) or fitted (as a discontinuity). 
The artificial viscosity term required to smooth out the discontinuity is 
usually of the same order as the mesh size (because of large, but finite, 
gradients of the solution in the shock region, even if higher order schemes are 
used). The same comment holds for finite elements. On the other hand, the 
discontinuous finite element approximation of shock waves proved to be efficient 
in nonlinear elasticity (see ref. 20). Here we will describe a finite element 
analogue for the shock fitting procedure used by Hafez and Cheng (ref. 21). 

Consider a shock element, as shown in figure 6. The Rankine-Hugoniot 
relations under the transonic small disturbance assumptions are given in 
equations (13) and (14). 

The first relation can be derived actually from the weak solution ad- 
mitted by TSDE, while the second is consistent with the irrotationality 
condition, which is equivalent to 0 . The equation for the nodal value 

at i - 1 will not be affected. The equation at i , however, will be 
different since only the contribution of segment II downstream of the shock 
will be considered. To the first order of accuracy, knowing and > 
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we know the condition upstream of the shock. can be determined according 

to relation (2) and in terms of 0 . . The rignthand side (dx/dy)2 may be 

evaluated from a previous iteration as the average of the slope of the shock 
in the adjacent elements. If this term is neglected, the scheme will reduce 
to the shock point operator, as discussed earlier. The compatibility relation 
(3’) is satisfied by using linear shape functions in upstream and downstream 
segments. Thus, (in finite difference calculations) the introduction of 
shock relations will not make the system matrix singular or disturb the con- 
vergence of iterations. 


As an alternative approach, instead of altering the nodal equation at the 
shock point to admit the jump in (|)^ between i and i-1, according to 
equation (1), we may use the divergence theorem to obtain an integral relation 
as a consexTvation of mass over the element. The element equation will read 


// 


V’ gdA 


g 



•nds = 0 
- 1/2 4 )^ 



( 22 ) 


Bilinear shape functions in I thru IV (fig* 6) may be used with a jump in 
across the shock. Similarly, the irrotationality condition (existence of 
potential) implies zero vorticity over each element and, by Stokes theorem, 
circulation, namely 


where 



(23) 


zero 


So, as an alternative approach, relations (2) and (3) are replaced by 
relations (22) and (23). 


REMARKS AND COMMENTS 


Mixed Variational Principles 


Note that higher order shape functions, namely. Hermits cubics, lead 
to equations (20) and (20’) for and at the nodes. The resulting 

algebraic equations can be considered as finite difference approximations of 
two differential equations: the first is the TSDE (1), and the second is the 
x-derivative of the TSDE. Instead, the problem can be formulated in terms of 
two unknown functions (fi and u , where c() is governed by the TSDE and u 
is governed by a comparability relation u = cj) . A mixed variational principle 
(in terms of (j) and u ), together with a dual iterative procedure for TSDE, 
is studied in a separate paper where the merits and the efficiency of the new 
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method is assessed. 


Weighted Residual Methods 


Chan and Brashears (ref. 9) used least squares to solve TSDE. Straight- 
forward application of the method fails (the solution diverges), so results 
can be obtained by changing the system matrix. The element matrices are 
constructed in the usual manner. Before assembling the element matrices into 
the system matrix, the rows corresponding to the nodes along the upstream side 
of any element in the supersonic zone are zeroed out; hence, no upstream effect 
is allowed there. Applying a similar procedure using the Galerkin method and 
cubic elements in the x-direction gives 


or 




H 


11 


(M -^) 

ax '■ ax'^ 


+ K 


1 w 


dx = 0 


(24) 


“ll “l2 


■{ 01 }' 



6i2 

”{ 



o 

4-1 

ttai a22 


{ 0 } 

L J 


_^21 

^22_ 

{ 

0 } 
XQ J 




(24») 


Note that equations (24) and (24’) differ from equations (20) and (20’) since 
different weighting functions are used. 


Type-Insensitive Methods 


In our method, a different strategy is adopted in subsonic and supersonic 
regions. A unified, type-insensitve method may be simpler, but not efficient, 
since different requirements in each region must be satisified simultaneously. 

To obtain such a procedure, the steady problem is embedded in a higher 
dimensional space, where the problem is more amenable for analysis. The extra 
dimension may have a physical meaning^ as in the unsteady (time-dependent) 
method or may be just a mathematical trick, like the use of -complex character- 
istics or any parameter as in the method of parametric differentiation. Also, 
extra dependent variables may be used, as in the mixed variational principle. 
The usefulness of these imbedding techniques depends on how fast the limit 
solution will be obtained. As an example of a unified procedure, consider the 
TSDE in the form of a system of first order equations, 


or 

C' .:) 


K^u = v = f = K-u 

£ X y £ 

u - v = g 

y X 

{:); C. “)(:), ■ 0 
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For cases where was a linear function of y , Friedrichs (ref. 22) 

and Chu (ref. 23) found a transformation that put this system into a positive 
symmetric form. As shown by Lesaint (ref. 24) and reported by Levanthal and 
Aziz (ref. 25), the finite element method can be applied successfully using 
this transformation. In general, however, such a transformation may not exist. 
Nevertheless, if the problem is considered as the asymptotic limit of an un- 

/f\ /"ct 0\ ful 

steady problem, where the vector Igl is replaced by Iq 3 / \v situation is 

different. The modified system' is symmetric and hyperbolic. Unlike the 
equilibrium equations, for symmetric hyperbolic equations positivity could always 
be attained by a simple transformation, as shown by Friedrichs (ref. 22). For 
such a modified system, no special treatment for subsonic and supersonic 
regions is needed. 

However, based on the finite difference calculations of the Euler 
equations, where centered differences are used everywhere in space, this 
"iterative” procedure may be slow. On the other hand, it seems that efficient 
applications of finite element methods to the full potential equation may 
require such imbedding techniques (artificial time-dependent and viscosity 
terms) , 


CONCLUSIONS 


Applications of a finite element approach to transonic flow problems have 
been discussed. Only small disturbance equations with streamlines almost 
parallel to the x-axis (hence, the nodes are located along x = constant lines 
in the supersonic region) have been considered. Currently, computations of a 
simple numerical example are underway. Extension of this approach to the full 
potential equation is possible as long as the direction of the flow in the 
supersonic region is almost known a priori. 
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Figure 1.- Mixed flow problems. 
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Figure 2.- Element for calculations using finite element 
in space, finite difference in time. 
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Figure 3.- Elements for calculations using finite element in space 

and time. 
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Figure 5.- Finite element assembly modified according to the inverse 

shooting method. 



Figure 6.- Shock element in finite element scheme. 
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